TP Integrador - Aprendizaje Estadístico

Enunciado

En el archivo sat.txt se tienen datos de valores de espectros de pixels en una imagen de satélite utilizados para la predicción de la clase de suelo. Realice un análisis de los datos utilizando diferentes métodos y compárelos mediante Validación Cruzada. Luego aplíquelos a la muestra de test sat.tst y analice los resultados obtenidos.

Descripción de los Datos

Propósito

La base de datos contiene los valores multi-espectrales de píxeles en regiones de 3x3 en una imágen satelital, y la clasificación asociada con el píxel central de cada región. El objetivo es predecir esta clasificación dados los valores multi-espectrales.

Detalle

Esta base de datos se generó a partir del escáner multi-espectral LANDSAT. Un marco de imágenes Landsat MSS consta de cuatro imágenes digitales de la misma escena en diferentes bandas espectrales. Dos de estos están en la región visible (correspondiente aproximadamente a las regiones verdes y rojas del espectro visible) y dos están en el infrarrojo. Cada píxel es una palabra binaria de 8 bits, con 0 correspondiente al negro y 255 a blanco. La resolución espacial de un píxel es de unos 80 m x 80m. Cada imagen contiene 2340 x 3380 de tales píxeles.

La base de datos es una subárea (pequeña) de una escena, que consta de 82 x 100 píxeles. Cada línea de datos corresponde a una región de píxeles cuadrados 3x3 completamente contenida dentro de la subárea 82x100. Cada línea contiene los valores de píxeles en las cuatro bandas espectrales (convertidas en ASCII) de cada uno de los 9 píxeles en la región 3x3 y un número que indica la etiqueta de clasificación del píxel central.

La clase de cada píxel se codifica como un número, siendo estos:

1: suelo rojo
2: cultivo de algodón
3: suelo gris
4: suelo gris húmdeo
5: suelo con rastrojo de vegetación
6: mezcla de suelos
7: suelo gris muy húmedo
- N.B: No hay ejemplos con la clase 6 en este conjunto de datos (fueron removidos debido a dudas respecto a la validez de esta clase). Los datos se dan en orden aleatorio y se han eliminado ciertas líneas de datos para que no pueda reconstruir la imagen original a partir de este conjunto de datos.

En cada línea de datos, los cuatro valores espectrales para el píxel superior izquierdo se dan primero, seguidos por los cuatro valores espectrales para el píxel medio superior y luego los del píxel superior derecho, y así sucesivamente de modo que los píxeles se leen en secuencia de izquierda a derecha y de arriba a abajo. Por lo tanto, los cuatro valores espectrales para los píxeles centrales están dados por los atributos 17,18,19 y 20. Si lo desea, puede usar solo estos cuatro atributos, ignorando los demás.

DATOS TOTALES → set de entrenamiento: 4435 set de prueba: 2000
CANTIDAD DE ATRIBUTOS → 36 (4 bandas espectrales x 9 píxeles en una región)
ATRIBUTOS → numéricos en el rango {0,255}
CLASES → 6 en total: {1,2,3,4,5,7}

Resolución

Obtención del Set de Datos

Dado que el enunciado menciona que pueden utilizarse solo los datos espectrales correspondientes al píxel central, buscamos extraerlos, junto con su clasificación, a una sola dataframe.

file_trn = "SAT_trn.txt"
data_trn <- read.table(file_trn, sep = "", header=F)
sat_trn_df <- data_trn      %>% 
  dplyr::select(17:20,37)  %>% 
  rename(B1=V17,
         B2=V18,
         B3=V19,
         B4=V20,
         C=V37)
file_tst = "SAT_tst.txt"
data_tst <- read.table(file_tst, sep = "", header=F)
sat_tst_df <- data_tst      %>% 
  dplyr::select(17:20,37)  %>% 
  rename(B1=V17,
         B2=V18,
         B3=V19,
         B4=V20,
         C=V37)

sat_trn_df$C[sat_trn_df$C == 7] <- 6
sat_tst_df$C[sat_tst_df$C == 7] <- 6
sat_trn_df_numeric <- data.frame(sat_trn_df)
sat_tst_df_numeric <- data.frame(sat_tst_df)
class_names <- c("rojo", "algodón", "gris","gris húmedo", "vegetación", "gris muy húmedo")
class_values <- c(1,2,3,4,5,6)
class_map <- setNames(class_names, class_values)
sat_trn_df$C <- class_map[sat_trn_df$C]
sat_tst_df$C <- class_map[sat_tst_df$C]

Visualización de los Datos

Contenido del Set de Datos

B1	B2	B3	B4	C
92	112	118	85	gris
84	103	104	81	gris
84	99	104	78	gris
84	99	104	81	gris
76	99	104	81	gris
76	99	108	85	gris
80	112	118	88	gris
80	107	113	85	gris
76	91	104	74	gris húmedo
76	95	100	78	gris húmedo

Podemos comprobar que tenemos valores acotados entre 0 y 255 para cuatro bandas espectrales y una categoría de suelo asignada en cada fila. Los contenidos pueden visualizarse mejor utilizando un gráfico de coordenadas paralelas.

Donde cada linea representa una observación para las 4 bandas. La normalización de los valores no es necesaria ya que \(B_1:B_4\) estan en la misma escala. Podemos ver que la distribución de suelo en los datos de entrenamiento y prueba es practicamente la misma.

Construcción de Modelos

sat_trn_df <- data.frame(sat_trn_df_numeric)
sat_tst_df <- data.frame(sat_tst_df_numeric)
sat_trn_df$C <- as.factor(sat_trn_df$C)
sat_tst_df$C <- as.factor(sat_tst_df$C)
set.seed(20)

Árbol de Decisión

Librería Tree

# Entrenamiento del Modelo
n <- nrow(sat_trn_df)
sat_tree <- tree(C ~ B1 + B2 + B3 + B4, 
                 data= sat_trn_df, 
                 control=tree.control(
                   nobs= nrow(sat_trn_df), 
                   mincut = 1, 
                   minsize = 2, 
                   mindev = 0)
                 )
# Resultados
sat_tree_summary <- summary(sat_tree)
EEnt <- sat_tree_summary$misclass[1] / n
sat_tree_summary

## 
## Classification tree:
## tree(formula = C ~ B1 + B2 + B3 + B4, data = sat_trn_df, control = tree.control(nobs = nrow(sat_trn_df), 
##     mincut = 1, minsize = 2, mindev = 0))
## Number of terminal nodes:  761 
## Residual mean deviance:  0.1723 = 633 / 3674 
## Misclassification error rate: 0.04014 = 178 / 4435

sat_tree_CV = cv.tree(sat_tree, FUN = prune.misclass)
size <- sat_tree_CV$size
dev <- sat_tree_CV$dev
sp <-ggplot(data=data.frame(size,dev), aes(x=size, y=dev)) +
  geom_line(linetype="solid", color="blue", size=1.2)+
  geom_point(color="red", size=3)+
  labs(x= "Tamaño del Árbol", y = "Error de CV",title="Poda de Complejidad de Costos")
sp + scale_x_continuous(trans='log10')

# Agrego código para obtener con R el mejor Tamaño
best_size <- which.min(sat_tree_CV$dev) 
best_size

## [1] 39

# PQ Harcodea el 13?
sat_tree_prune = prune.misclass(sat_tree, best = 13)
plot(sat_tree_prune)
text(sat_tree_prune, pretty = 0)

CVErr_tree <- sat_tree_CV$dev[39] / n
CVErr_tree

## [1] 0.1815107

Librería RPart

set.seed(20)
# Entrenamiento del Modelo con 10-fold CV  
sat_rpart <-  rpart(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method = "class", cp=-1,xval=10)
sat_rpart_summary <- printcp(sat_rpart)

## 
## Classification tree:
## rpart(formula = C ~ B1 + B2 + B3 + B4, data = sat_trn_df, method = "class", 
##     cp = -1, xval = 10)
## 
## Variables actually used in tree construction:
## [1] B1 B2 B3 B4
## 
## Root node error: 3363/4435 = 0.75829
## 
## n= 4435 
## 
##             CP nsplit rel error  xerror      xstd
## 1   2.6167e-01      0   1.00000 1.00000 0.0084779
## 2   2.5335e-01      1   0.73833 0.76896 0.0097636
## 3   1.2905e-01      2   0.48498 0.48498 0.0095487
## 4   6.3634e-02      3   0.35593 0.35742 0.0088020
## 5   1.7544e-02      4   0.29230 0.29379 0.0082400
## 6   1.4570e-02      5   0.27475 0.28189 0.0081181
## 7   6.6905e-03      7   0.24561 0.24234 0.0076694
## 8   5.9471e-03     10   0.22034 0.23253 0.0075467
## 9   5.6497e-03     11   0.21439 0.22629 0.0074660
## 10  4.7577e-03     13   0.20309 0.21885 0.0073673
## 11  4.6090e-03     14   0.19833 0.20934 0.0072364
## 12  3.2709e-03     16   0.18912 0.20636 0.0071945
## 13  1.6354e-03     17   0.18585 0.19774 0.0070698
## 14  1.3381e-03     19   0.18258 0.19625 0.0070479
## 15  1.1894e-03     21   0.17990 0.19536 0.0070346
## 16  1.0903e-03     22   0.17871 0.19477 0.0070258
## 17  8.9206e-04     25   0.17544 0.19715 0.0070610
## 18  6.9382e-04     31   0.17009 0.19715 0.0070610
## 19  5.9471e-04     34   0.16800 0.19715 0.0070610
## 20  4.9559e-04     48   0.15968 0.19685 0.0070567
## 21  4.4603e-04     51   0.15819 0.19387 0.0070124
## 22  3.9647e-04     55   0.15641 0.19328 0.0070035
## 23  2.9735e-04     58   0.15522 0.19328 0.0070035
## 24  1.4868e-04     61   0.15433 0.20012 0.0071047
## 25  8.4958e-05     63   0.15403 0.20101 0.0071176
## 26  7.4338e-05     70   0.15343 0.20071 0.0071133
## 27  0.0000e+00     74   0.15314 0.19952 0.0070960
## 28 -1.0000e+00    197   0.15314 0.19952 0.0070960

plotcp(sat_rpart)

which.min(sat_rpart$cptable[,4])

## 22 
## 22

cp<-sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"CP"]
cp

## [1] 0.0003964714

CVErr_rpart = sat_rpart$frame[1,'dev'] / n * sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"xerror"]
CVErr_rpart

## [1] 0.1465614

poda_rpart<-prune(sat_rpart,cp=cp)
rpart.plot(poda_rpart, cex = 0.5, extra = 0)

Random Forest

set.seed(20)
sat_rf<- randomForest(C~., data= sat_trn_df, 
                     mtry=sqrt(ncol(sat_trn_df)-1), importance=TRUE)

importance(sat_rf)

##            1        2         3        4         5        6
## B1 227.87903 29.96135 362.16076 58.21353  50.21856 48.37936
## B2 140.86987 43.00905  53.51638 37.62347 112.03571 41.93234
## B3  30.04423 10.40657  24.77468 15.11574  28.56435 32.31841
## B4  57.76868 24.80617  34.26451 35.82601  52.16778 74.12139
##    MeanDecreaseAccuracy MeanDecreaseGini
## B1            186.14425        1186.7218
## B2            119.11363         914.1809
## B3             48.88156         425.2869
## B4             91.94379         853.1770

varImpPlot(sat_rf)

CVErr_rf <- unname(sat_rf$err.rate[nrow(sat_rf$err.rate),1])
CVErr_rf

## [1] 0.1492672

Vecinos Más Cercanos

set.seed(20)
sat_knn <- train(C ~ ., data = sat_trn_df, 
                   method = "knn",
                   preProcess = c("center","scale"),
                   trControl =  trainControl(method = "cv", number = 10),
                   tuneGrid = expand.grid(k = 1:100),
                   metric = "Accuracy")
results <- sat_knn$results
ggplot(results, aes(x = k, y = Accuracy)) + 
  geom_line(color = "steelblue", size = 1.2) + 
  xlab("Valor de K") + 
  ylab("Exactitud") + 
  ggtitle("Rendimiento para diferentes valores de K") +
  theme_bw() + 
  theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
        axis.title.x = element_text(face = "bold", size = 12),
        axis.title.y = element_text(face = "bold", size = 12),
        axis.text = element_text(size = 10),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

best_k <- results$k[which.max(results$Accuracy)]
best_k

## [1] 43

CVErr_KNN <- 1-unname(sat_knn$results[best_k,"Accuracy"])
CVErr_KNN

## [1] 0.1400425

Regresión Logística Multinomial

summary(sat_mlr)

## Call:
## nnet::multinom(formula = .outcome ~ ., data = dat, decay = param$decay)
## 
## Coefficients:
##   (Intercept)        B1          B2          B3         B4
## 2   -9.921154 0.5867093 -0.63317495  0.29943995 -0.1122932
## 3  -36.648978 0.7617752 -0.02673514  0.02795903 -0.2534971
## 4  -10.253548 0.5886442 -0.07999703 -0.02004337 -0.2854203
## 5    7.594031 0.4558313 -0.52173889  0.21832973 -0.2099831
## 6    4.599893 0.5678520 -0.18940106 -0.05394639 -0.2887664
## 
## Std. Errors:
##   (Intercept)         B1         B2         B3         B4
## 2   1.6977543 0.05040409 0.04477873 0.05284647 0.04508065
## 3   1.2832527 0.03391253 0.03743092 0.04631104 0.04434272
## 4   0.8182467 0.03437754 0.03467115 0.04409015 0.04040293
## 5   0.9108139 0.03161938 0.03489501 0.04077884 0.03601165
## 6   0.7841282 0.03491517 0.03359844 0.04276092 0.03778148
## 
## Residual Deviance: 3755.425 
## AIC: 3805.425

#sat_mlr <- nnet::multinom(C ~ ., data = sat_trn_df)
#summary(sat_rl)
# Predict the class labels for the training data
confusion_matrix <- confusionMatrix(predict(sat_mlr, newdata = sat_trn_df), sat_trn_df$C)

CVErr_mlr2 <- 1 - unname(confusion_matrix$overall['Accuracy'])
CVErr_mlr2

## [1] 0.1569335

predicted_mlr <- predict(sat_mlr, newdata = sat_trn_df)

CVErr_mlr <- mean(predicted_mlr != sat_trn_df$C)
CVErr_mlr

## [1] 0.1569335

Análisis Discriminante Lineal

set.seed(20)
sat_LDA <- stepclass(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method="lda", criterion = "CR", direction = "both", improvement <- 0.01, fold = 10)

##  `stepwise classification', using 10-fold cross-validated correctness rate of method lda'.

## 4435 observations of 4 variables in 6 classes; direction: both

## stop criterion: improvement less than 1%.

## correctness rate: 0.56054;  in: "B2";  variables (1): B2 
## correctness rate: 0.79776;  in: "B1";  variables (2): B2, B1 
## correctness rate: 0.82616;  in: "B3";  variables (3): B2, B1, B3 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##        0.00        0.00        3.13

sat_LDA

## method      : lda 
## final model : C ~ B1 + B2 + B3
## <environment: 0x000001f8ba1e8210>
## 
## correctness rate = 0.8262

CVErr_LDA <- 1-unname(sat_LDA$result.pm["crossval.rate"])

Análisis Discriminante Cuadrático

set.seed(20)
sat_QDA <- stepclass(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method="qda", criterion = "CR", direction = "both", improvement <- 0.01, fold = 10)

##  `stepwise classification', using 10-fold cross-validated correctness rate of method qda'.

## 4435 observations of 4 variables in 6 classes; direction: both

## stop criterion: improvement less than 1%.

## correctness rate: 0.58107;  in: "B4";  variables (1): B4 
## correctness rate: 0.80496;  in: "B1";  variables (2): B4, B1 
## correctness rate: 0.84849;  in: "B2";  variables (3): B4, B1, B2 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##        0.00        0.00        3.44

sat_QDA

## method      : qda 
## final model : C ~ B1 + B2 + B4
## <environment: 0x000001f8bde91e90>
## 
## correctness rate = 0.8485

CVErr_QDA <- 1-unname(sat_QDA$result.pm["crossval.rate"])

Comparación del Error de 10-fold CV

models <- c("tree", "rpart", "RF", "KNN", "RLM","LDA","QDA")
CVErrs <- c(CVErr_tree,CVErr_rpart,CVErr_rf,CVErr_KNN,CVErr_mlr,CVErr_LDA,CVErr_QDA)

CVErr_df <- data.frame(Label = models, Value = CVErrs)

ggplot(CVErr_df, aes(x = Label, y = Value)) + 
  geom_bar(stat = "identity", color = "black", fill = "red") + labs(x = "Modelos", y = "Error de Validación Cruzada K=10", title = "Comparación de Precisión de Modelos") +
  theme(plot.title = element_text(hjust = 0.5))

#conf_mat <- caret::confusionMatrix(pred_nnet, sat_tst_df$C)
#conf_mat

#conf_matrix(sat_tst_df$C,pred_nnet)

#probs <- predict(model_nnet, newdata = sat_tst_df, type = "probs")
#ROC <- multiclass.roc(response = sat_tst_df$C, predictor = probs)

#auc(ROC)

Evaluación de los Modelos contra el Set de Entrenamiento

Árbol de Clasificación

Librería Tree

Librería Rpart

Random Forest

Regresión Logística Multinomial

Vecinos Más Cercanos

Análisis Discriminante Lineal

Análisis Discriminante Cuadrático

ANEXO: Análisis Para el Set de Datos Completo

Queremos comprobar:

Que el píxel central es en efecto el principal predictor de la clase de la celda
Que despreciar los píxeles adyacentes al píxel central no produce una reducción significativa en el error de clasificación.

sat_trn_df_full <- data_trn
sat_trn_df_full$V37 <- as.factor(sat_trn_df_full$V37)

n <- nrow(sat_trn_df_full)
full_tree <- tree(V37 ~., data= sat_trn_df_full, control=tree.control(nobs= nrow(sat_trn_df_full), mincut = 1, minsize = 2))
full_tree_summary <- summary(full_tree)
EEnt <- full_tree_summary$misclass[1] / n
#full_tree_summary

set.seed(20)
full_tree_CV = cv.tree(full_tree, FUN = prune.misclass)
full_size <- full_tree_CV$size
full_dev <- full_tree_CV$dev
#ggplot(data=data.frame(full_size,full_dev), aes(x=full_size, y=full_dev)) +
#  geom_line(linetype="solid", color="blue", size=1.2)+
#  geom_point(color="red", size=3)+
#  labs(x= "Tamaño del Árbol", y = "Error de CV",title="Poda de Complejidad de Costos")

prune_full = prune.misclass(full_tree, best = 10)
plot(prune_full)
text(prune_full, pretty = 0)

set.seed(20)
full_tree_2 <-  rpart(V37 ~., data= sat_trn_df_full, method = "class", cp=-1,xval=10)
full_tree_2_summary <- printcp(full_tree_2)

## 
## Classification tree:
## rpart(formula = V37 ~ ., data = sat_trn_df_full, method = "class", 
##     cp = -1, xval = 10)
## 
## Variables actually used in tree construction:
##  [1] V1  V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V2  V20 V21 V22 V23 V24 V25 V26
## [20] V27 V28 V29 V3  V30 V31 V32 V33 V34 V35 V36 V4  V5  V6  V7  V8  V9 
## 
## Root node error: 3363/4435 = 0.75829
## 
## n= 4435 
## 
##             CP nsplit rel error  xerror      xstd
## 1   0.26167113      0   1.00000 1.00000 0.0084779
## 2   0.25334523      1   0.73833 0.76896 0.0097636
## 3   0.12905144      2   0.48498 0.48498 0.0095487
## 4   0.06363366      3   0.35593 0.35742 0.0088020
## 5   0.01754386      4   0.29230 0.29497 0.0082518
## 6   0.01486768      5   0.27475 0.28338 0.0081337
## 7   0.01159679      6   0.25989 0.25453 0.0078153
## 8   0.00683913      7   0.24829 0.24591 0.0077128
## 9   0.00669045      9   0.23461 0.23877 0.0076253
## 10  0.00624442     13   0.20458 0.23699 0.0076030
## 11  0.00490633     15   0.19209 0.21588 0.0073269
## 12  0.00356824     17   0.18228 0.21053 0.0072531
## 13  0.00297354     18   0.17871 0.20607 0.0071903
## 14  0.00267618     19   0.17574 0.20517 0.0071776
## 15  0.00237883     21   0.17038 0.20369 0.0071563
## 16  0.00223015     22   0.16800 0.19923 0.0070916
## 17  0.00208147     24   0.16354 0.19923 0.0070916
## 18  0.00178412     26   0.15938 0.19744 0.0070654
## 19  0.00168500     27   0.15760 0.19685 0.0070567
## 20  0.00163544     30   0.15254 0.19685 0.0070567
## 21  0.00148677     32   0.14927 0.19506 0.0070302
## 22  0.00133809     36   0.14332 0.19387 0.0070124
## 23  0.00118941     40   0.13797 0.19447 0.0070213
## 24  0.00104074     50   0.12608 0.19387 0.0070124
## 25  0.00089206     55   0.12073 0.19120 0.0069721
## 26  0.00077312     57   0.11894 0.19060 0.0069631
## 27  0.00074338     62   0.11508 0.19090 0.0069676
## 28  0.00059471     64   0.11359 0.19090 0.0069676
## 29  0.00049559     66   0.11240 0.19179 0.0069811
## 30  0.00029735     69   0.11091 0.19269 0.0069946
## 31  0.00019824     74   0.10943 0.19804 0.0070742
## 32  0.00014868     77   0.10883 0.19893 0.0070873
## 33  0.00000000     83   0.10794 0.19982 0.0071003
## 34 -1.00000000    157   0.10794 0.19982 0.0071003

full_tree_2_NP_pred <- predict(full_tree_2, sat_trn_df_full, type = "class")
# Poda
#full_tree_2$cptable 
#plotcp(full_tree_2) #grafico
#which.min(full_tree_2$cptable[,4])
cp<-full_tree_2$cptable[which.min(full_tree_2$cptable[,"xerror"]),"CP"]
poda<-prune(full_tree_2,cp=cp)
rpart.plot(poda, cex = 0.5, extra = 0)

set.seed(20)
rf_full<- randomForest(V37~., data= sat_trn_df_full, 
                     mtry=sqrt(ncol(sat_trn_df_full)-1), importance=TRUE)

#importance(rf_full)
varImpPlot(rf_full)

unname(rf_full$err.rate[nrow(rf_full$err.rate),1])

## [1] 0.08568207

TP Integrador - Aprendizaje Estadístico

Santiago Bertero - Lautaro De Lucia

2023-01-17

Enunciado

Descripción de los Datos

Propósito

Detalle

Resolución

Obtención del Set de Datos

Visualización de los Datos

Contenido del Set de Datos

Construcción de Modelos

Árbol de Decisión

Librería Tree

Librería RPart

Random Forest

Vecinos Más Cercanos

Regresión Logística Multinomial

Análisis Discriminante Lineal

Análisis Discriminante Cuadrático

Comparación del Error de 10-fold CV

Evaluación de los Modelos contra el Set de Entrenamiento

Árbol de Clasificación

Librería Tree

Librería Rpart

Random Forest

Regresión Logística Multinomial

Vecinos Más Cercanos

Análisis Discriminante Lineal

Análisis Discriminante Cuadrático

ANEXO: Análisis Para el Set de Datos Completo